function mpr = getMPR(output,theta2,AVarTheta,nx,h0RegimeOn,hxRegimeOn)

% Version 1: MPR = (lambda0 + lambdax)
% Version 2: MPR = (lambda0 + lambdax)*inv(sigma)
version = 2;    

[h0_1,h0_2,hx_1,hx_2,sigma] = unfoldTheta2_threshold_macro(theta2,nx);
if version == 1
    lambda0_1 = sigma\(h0_1-output.model.phi*output.model.mu);
    lambda0_2 = sigma\(h0_2-output.model.phi*output.model.mu);
    lambdax_1 = sigma\(hx_1-(eye(nx)-output.model.phi));
    lambdax_2 = sigma\(hx_2-(eye(nx)-output.model.phi));
elseif version == 2
    lambda0_1 = h0_1-output.model.phi*output.model.mu;
    lambda0_2 = h0_2-output.model.phi*output.model.mu;
    lambdax_1 = (hx_1-(eye(nx)-output.model.phi));
    lambdax_2 = (hx_2-(eye(nx)-output.model.phi));
end

if version == 1
    % Computing the standard errors by simulation (ignoring uncertainty about
    % theta1, which is ok for T going to infinitey because then N*T goes much
    % faster to infinity given N also goes to infinity)
    numSim = 50000;
    rng(1,'twister');
    lambda0_1Draw = zeros(nx,numSim);
    lambda0_2Draw = zeros(nx,numSim);
    lambdax_1Draw = zeros(nx,nx,numSim);
    lambdax_2Draw = zeros(nx,nx,numSim);
    namesTheta2Reduced = AVarTheta.names;
    numTheta2Reduced   = length(namesTheta2Reduced);
    for i=1:numTheta2Reduced
        theta2Reduced.(namesTheta2Reduced{i}) = theta2.(namesTheta2Reduced{i});
    end
    [Stheta2Reduced,check] = chol(AVarTheta.VarRobust,'lower');
    if check ~= 0
        mpr.lambda0_1       = lambda0_1;
        mpr.lambda0_2       = lambda0_2;
        mpr.lambdax_1       = lambdax_1;
        mpr.lambdax_2       = lambdax_2;
        mpr.SE = NaN;
        return
    end
    for s=1:numSim
        theta2ReducedDrawValues = struct2array(theta2Reduced)'+Stheta2Reduced*randn(numTheta2Reduced,1);
        theta2ReducedDraw       = values2struct(theta2ReducedDrawValues,namesTheta2Reduced);
        if h0RegimeOn == 0
            for i=1:nx
                theta2ReducedDraw.(['h0',num2str(i),'_2']) = theta2ReducedDraw.(['h0',num2str(i),'_1']);
            end
        end
        if hxRegimeOn == 0
            for i=1:nx
                for j=1:nx
                    theta2ReducedDraw.(['hx',num2str(i),num2str(j),'_2']) = theta2ReducedDraw.(['hx',num2str(i),num2str(j),'_1']);
                end
            end
        end
        [h0_1,h0_2,hx_1,hx_2,sigma] = unfoldTheta2_threshold_macro(theta2ReducedDraw,nx);
        % Intercepts for MPR
        lambda0_1Draw(:,s) = sigma\(h0_1-output.model.phi*output.model.mu);
        lambda0_2Draw(:,s) = sigma\(h0_2-output.model.phi*output.model.mu);
        %lambda0_1Draw(:,s) = h0_1-output.model.phi*output.model.mu;
        %lambda0_2Draw(:,s) = h0_2-output.model.phi*output.model.mu;
        % Loadings for MPR
        lambdax_1Draw(:,:,s) = sigma\(hx_1-(eye(nx)-output.model.phi));
        lambdax_2Draw(:,:,s) = sigma\(hx_2-(eye(nx)-output.model.phi));
        %lambdax_1Draw(:,:,s) = (hx_1-(eye(nx)-output.model.phi));
        %lambdax_2Draw(:,:,s) = (hx_2-(eye(nx)-output.model.phi));
    end
    lambda0_1_SE = std(lambda0_1Draw,0,2);
    lambda0_2_SE = std(lambda0_2Draw,0,2);
    lambdax_1_SE = std(lambdax_1Draw,0,3);
    lambdax_2_SE = std(lambdax_2Draw,0,3);
elseif version == 2
    theta2se = struct();
    for i=1:length(AVarTheta.names)
        name = AVarTheta.names{i};
        theta2se.(name) = sqrt(AVarTheta.VarRobust(i,i));
    end
    [h0_1se,h0_2se,hx_1se,hx_2se] = unfoldTheta2_threshold_macro(theta2se,nx);
    lambda0_1_SE = h0_1se;
    lambda0_2_SE = h0_2se;
    lambdax_1_SE = hx_1se;
    lambdax_2_SE = hx_2se;
end
% Computing t-stats
lambda0_1_tstat = lambda0_1./lambda0_1_SE;
lambda0_2_tstat = lambda0_2./lambda0_2_SE;
lambdax_1_tstat = lambdax_1./lambdax_1_SE;
lambdax_2_tstat = lambdax_2./lambdax_2_SE;

%% saving results
% nameLambda0_1 = {};
% nameLambda0_2 = {};
% for i=1:nx
%     nameLambda0_1 =  [nameLambda0_1 {['lambda0',num2str(i),'_1']}];
%     nameLambda0_2 =  [nameLambda0_2 {['lambda0',num2str(i),'_2']}];
% end
% nameLambdax_1 = {};
% nameLambdax_2 = {};
% for i=1:nx
%     for j=1:nx
%         nameLambdax_1 =  [nameLambdax_1 {['lambdax',num2str(i),num2str(j),'_1']}];
%         nameLambdax_2 =  [nameLambdax_2 {['lambdax',num2str(i),num2str(j),'_2']}];
%     end
% end
mpr.lambda0_1       = lambda0_1;
mpr.lambda0_2       = lambda0_2;
mpr.lambdax_1       = lambdax_1;
mpr.lambdax_2       = lambdax_2;
mpr.lambda0_1_SE    = lambda0_1_SE;
mpr.lambda0_2_SE    = lambda0_2_SE;
mpr.lambdax_1_SE    = lambdax_1_SE;
mpr.lambdax_2_SE    = lambdax_2_SE;
mpr.lambda0_1_tstat = lambda0_1_tstat;
mpr.lambda0_2_tstat = lambda0_2_tstat;
mpr.lambdax_1_tstat = lambdax_1_tstat;
mpr.lambdax_2_tstat = lambdax_2_tstat;
end